#include "FEMSpace.h"
#include <typeinfo>
#include <math.h>
#include <ctime>
#include "Error.h"
#define pi 4.0*atan(1.0)
double f(double *p)
{
    return 2 * pi*  pi* sin(pi * p[0]) * sin(pi * p[1]);
}

double bc(double *p)
{
    return sin(pi * p[0]) * sin(pi * p[1]);
}

double f2(double *p)
{
    return 2 * pi*  pi* cos(pi * p[0]) * cos(pi * p[1]);
}

double bc2(double *p)
{
    return cos(pi * p[0]) * cos(pi * p[1]);
}
int main()
{
    //g++ -o main test_possion.cpp -std=c++11 -I /usr/include/eigen3/
    RectangleDomain* r = new RectangleDomain({{0,0},{1,0},{1,1},{0,1}});
    std::vector<Boundary<2> > B = r->boundary();
    int segmentx = 16;
    int segmenty = segmentx;
    Mesh<2>* m = new Q1Mesh(r,{segmentx,segmenty});
    Element<2>* e = new Quadrilateral_1_Element();
    Equation<2>* equ = new PossionEquation<2>();
    equ ->SetRightHandsTermFunction(f2);
    BoundaryFunction<2> * bf = new Dirichlet<2>(bc2,B);
    BoundaryCondition<2> bfc;
    bfc.add(bf);
    Possion_2D possionproblem(m,e,equ,bfc);
    possionproblem.AssembleStiffMatrix();
    possionproblem.AssembleRightHandsTerm();
    possionproblem.DealWithBoubdaryCondition();
    possionproblem.Solve();
}
